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Abstract. In the framework of the canonical model of hy- 
drodynamics, where fluid is assumed to be ideal and incom- 
pressible, waves are potential, two-dimensional, and sym- 
metric, the authors have recently reported the existence of a 
new type of gravity waves on deep water besides well studied 
Stokes waves (Lukomsky et al., 2002b). The distinctive fea- 
ture of these waves is that horizontal water velocities in the 
wave crests exceed the speed of the crests themselves. Such 
waves were found to describe irregular flows with stagnation 
point inside the flow domain and discontinuous streamlines 
near the wave crests. Irregular flows produce a simple model 
for describing the initial stage of the formation of spilling 
breakers when a localized jet is formed at the crest following 
by generating whitecaps. 

In the present work, a new highly efficient method for 
computing steady potential gravity waves on deep water 
is proposed to examine the above results in more detail. 
The method is based on the truncated fractional approxima- 
tions for the velocity potential in terms of the basis func- 
tions l/(l — exp(y — V — ix)) , yo being a free param- 
eter. The non-linear transformation of the horizontal scale 
x = x — isinx, < 7 < 1 , is additionally applied to 
concentrate a numerical emphasis on the crest region of a 
wave for accelerating the convergence of the series. Frac- 
tional approximations were employed for calculating both 
steep Stokes waves and irregular flows. For lesser computa- 
tional time, the advantage in accuracy over ordinary Fourier 
expansions in terms the basis functions cxp(n(y + ix)) was 
found to be from one to ten decimal orders depending on the 
wave steepness and flow parameters. 



1 Introduction 

From old times the wave motion of the ocean bewitched and 
extremely attracted the attention of mankind. Up to now the 
problem of understanding specific features of water waves 
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and their modelling represent a real challenge both from sci- 
entific and engineering points of view. Occurrence of ex- 
tremely large and steep ocean breaking waves imposes a haz- 
ard to fishing boats, ships, and off-shore oil facilities. To 
understand physical mechanisms that give rise to extreme 
breaking waves and to model them correctly it is necessary to 
gain detailed knowledge of the form and dynamics of steep 
water waves. 

The canonical problem about the propagation of surface 
waves on deep water (see Section 2) was the first essentially 
non-linear problem in hydrodynamics. Its analysis during 
almost two hundred years gave the origin to many fields of 
non-linear dynamics such as solitary waves, modulation in- 
stabilities, strange attractors, etc. Stokes (1847) was the first 
who considered surface waves of finite amplitude (Stokes 
waves). Stokes (1880) conjectured that such waves must 
have a maximal amplitude (the limiting wave) and showed 
the flow in this wave to be singular at the crest forming a 120° 
corner (the Stokes corner flow). Much later, Grant (1973) 
suggested that this singularity, for a wave that has not at- 
tained the limiting form, is located above the wave crest and 
forms a stagnation point with streamlines meeting at right 
angles. Longuet-Higgins and Fox (1978) proved this numer- 
ically after extending Stokes flows analytically outside the 
domain filled by fluid. The following question resulted: why 
the flow in the limiting Stokes wave has the 120° singularity 
instead of the 90° one, as in any wave with lesser ampli- 
tude? Because of this Grant (1973) conjectured that a con- 
tinuous approach to the limiting amplitude is possible only 
if the Stokes corner flow has several coalescing singularities. 
However, it has not yet clear where these multiple singulari- 
ties arise from. 

A new era in developing the theory of steep gravity waves 
started from the work of Longuet-Higgins (1975), where he 
found that many characteristics of gravity waves, such as 
speed, energy, and momentum, are not monotonic functions 
of the wave amplitude, as was assumed from Stokes, but at- 
tain total maxima and then drop before the limiting wave 
is reached. Longuet-Higgins and Fox (1977) constructed 
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asymptotic expansions for waves close to the 120°-cusped 
wave (almost highest waves) and showed that these depen- 
dences oscillate infinitely as the limiting wave is approached. 
Nevertheless, strict numerical verification of such oscilla- 
tions seems to be a real challenge up to the present time, 
with only the first relative maximum and minimum having 
been thoroughly investigated (see, e.g., Longuet-Higgins and 
Tanaka, 1997). 

Tanaka (1983) showed that gravity waves steeper than the 
wave with maximal total energy become unstable with re- 
spect to two-dimensional disturbances having the same pe- 
riod as an undisturbed wave (superharmonic instability). Jil- 
lians (1989) investigated the form of such instabilities and 
showed that they lead to wave overturning and breaking. 
The conjecture was made that wave breaking is a purely lo- 
cal phenomenon around the wave crest which, in the case 
of spilling breakers and more gently plunging breakers, oc- 
curs independently of the flow in the rest of a wave. Pro- 
ceeding with this idea Longuet-Higgins and Cleaver (1994) 
and Longuet-Higgins et al. (1994) suggested that superhar- 
monic instability results in the crests of almost highest Stokes 
waves to be unstable (crest instability). Longuet-Higgins and 
Tanaka (1997) strongly supported the conclusion that super- 
harmonic instabilities of Stokes waves are indeed crest insta- 
bilities. Finally, Longuet-Higgins and Dommermuth (1997) 
showed that crest instabilities lead (i) to wave overturning 
and breaking or (ii) to a smooth transition of a wave to a 
lower progressive wave having nearly the same total energy, 
followed by a return to a wave of almost the initial wave 
height. The latter fact generated a new question: what is the 
nature of such a transient phenomenon? A possible explana- 
tion would be found if superharmonic instability resulted in 
a bifurcation to a new solution, as usually takes place in non- 
linear dynamics. However, up to this time it was assumed 
that the Stokes solution is unique and free of bifurcations in 
keeping with the uniqueness argument of Garabedian (1965). 
The only bifurcation known to occur is the trivial one of a 
pure phase shift at the point of energy maximum (Tanaka, 
1985). 

The above results are all related to Stokes waves, for which 
the speed of fluid particles at the wave crests is smaller than 
the wave phase speed, equality being achieved for the limit- 
ing wave only. Thus, the traditional criterion for wave break- 
ing is that horizontal water velocities in the crest must ex- 
ceed the speed of the crest (Banner and Peregrine, 1993). 
Lukomsky et al. (2002a,b) have recently provided evidence 
(although numerical and not completely rigorous) for the ex- 
istence of a new family of two-dimensional irrotational sym- 
metric periodic gravity waves that satisfy the criterion of 
breaking. A stagnation point in the flow field of these waves 
is inside the flow domain, in contrast to the Stokes waves of 
the same wavelength. This makes streamlines exhibit discon- 
tinuity in the vicinity of the wave crests, with near-surface 
particles being jetted out from the flow. Because of this such 
waves and flows were called irregular (in contrast to regu- 
lar Stokes flows). Existence of irregular waves made the au- 
thors propose the following conjectures (see Lukomsky et al., 



2002b): (i) a 120° corner at the crest of the limiting Stokes 
wave seems to be formed by merging the stagnation points 
of a regular Stokes flow and an irregular flow; (ii) irregular 
flows are apparently the result of developing superharmonic 
or crest instabilities of regular Stokes flows. 

To calculate irregular flows Lukomsky et al. (2002a,b) 
used truncated Fourier expansions for the velocity poten- 
tial and the elevation of a free surface in the plane of 
spatial variables (a physical plane). Debiane and Kharif 
(see Gandzha et al., 2002) confirmed the existence of ir- 
regular waves using inverse plane Longuet-Higgins method 
(Longuet-Higgins, 1986), where the spatial coordinates are 
represented as Fourier series in velocity potential and stream 
function, the corresponding coefficients being evaluated by 
solving quadratic relations between them. Finally, Clam- 
ond (2003) also obtained irregular flows by applying his new 
renormalized cnoidal wave (RCW) approximation. Note, 
however, that at present these are only numerical and approx- 
imate proofs. 

Ordinary Fourier expansions used by Lukomsky et al. 
(2002a) become not efficient enough for approximating ir- 
regular waves and even Stokes waves close to the limiting 
one due to slow descending of the Fourier coefficients. Thus, 
a more efficient method is necessary for these tasks. Up 
to now the most precise and efficient way for calculating 
the properties of two-dimensional surface waves is Tanaka's 
method of the inverse plane (see Tanaka, 1983, 1986). The 
key idea of his method is to map the inverse plane into a 
unit circle by means of the Nekrasov transformation. Then 
boundary conditions are transformed to an integral equation, 
which is solved iteratively. The accuracy of obtained solu- 
tions is drastically improved by concentrating a numerical 
emphasis on the crest region using further transformation of 
variables. As a result, Tanaka's method is the only one be- 
ing capable of evaluating the second maximum of the phase 
speed and even further higher order extremums. In spite 
of all the advantages of Tanaka's method and his program, 
where it is implemented, we are interested in improving the 
methods of the physical plane since they can be applied for 
calculating 3D waves as well and can be generalized to the 
case of non-ideal and compressible fluid, in contrast to all the 
inverse plane methods. 

Thus, the purpose of this paper is to present a new method 
in the physical plane for calculating two-dimensional poten- 
tial steady progressive surface waves on the fluid of infinite 
depth (see Section 3). The method is based on the fractional 
Fourier approximation for the velocity potential recently in- 
troduced by the authors (Gandzha et al., 2002; Lukomsky 
et al., 2002c) and the non-linear transformation of the hori- 
zontal scale for concentrating a numerical emphasis near the 
wave crest. The first term of such a fractional Fourier ap- 
proximation was independently derived by Clamond (2003) 
and was called a renormalized cnoidal wave approximation. 

In Section 4, fractional Fourier approximations are applied 
for calculating regular and irregular flows. In Section 4.1, 
great numerical advantage of fractional approximations over 
ordinary Fourier approximations is demonstrated when cal- 
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Fig. 1. The laboratory frame of reference. 



culating almost highest Stokes waves. Although the accu- 
racy of the results is still less than the ones obtained us- 
ing Tanaka's method, proposed fractional Fourier approxi- 
mations have a potential to become almost as effective as 
the method of Tanaka. In Section 4.2, fractional approxima- 
tions allowed investigating the character of singularity of ir- 
regular flows in more detail. The profiles of irregular waves 
are shown to reveal the Gibbs phenomenon usually taking 
place when a discontinuous function or a continuous function 
with discontinuous derivatives are approximated by continu- 
ous truncated Fourier series. Moreover, even regular Stokes 
waves very close to the Stokes corner flow are also demon- 
strated to exhibit the similar Gibbs phenomenon in accor- 
dance with the observation of Chandler and Graham (1993). 
Concluding remarks are given in Section 5. 



2 The canonical model 

Consider the dynamics of potential two-dimensional periodic 
waves on the irrotational, inviscid, incompressible fluid with 
unknown free surface under the influence of gravity. Waves 
are assumed to propagate without changing their form from 
left to right along the x-axis with constant speed c relative 
to the motionless fluid at infinite depth (see Fig. 1). Gravity 
waves and related fluid flows are governed by the following 
set of equations 



§69 + ®yy = 0, 

(c - $g f + $1 + 2r] = c 2 , 
(c- $ fl )rye + $y = 0, 
$ = 0, $ y = 0, 



-oo < y < r](8); 

y = v(8); 



n(6) 

— oo 



(1) 

(2) 
(3) 
(4) 



where 9 = x — ct is the wave phase, y) is the velocity 
potential (the velocity is equal to V<P), 77(61) is the elevation 
of the unknown free surface, and y is the upward vertical 
axis such that y = is the still water level. Herein (1) is the 
Laplace equation in the flow domain, (2) is the dynamical 
boundary condition (the Bernoulli equation), (3) is the kine- 
matic boundary condition (no fluid crosses the surface), (4) 
is the condition that fluid is motionless at infinite depth. The 
dimensionless variables are chosen such that length and time 
are normalized by the wavenumber k and the frequency ^fgk 
of a linear wave, respectively, g being the acceleration due to 
gravity. In this case, the dimensionless wavelength A = 27r. 



Once the velocity potential and the wave phase speed are 
known, particle trajectories (streamlines) are found from the 
following differential equations: 



% y) ; 



(5) 



in the frame of reference moving together with the wave 
(the wave related frame of reference). In this frame of ref- 
erence, each streamline is characterized by a constant value 
of a stream function ip(6, y). The velocity potential and the 
stream function ^(9, y) = ip(9, y) + cy in the laboratory 
frame of reference are connected by means of the Cauchy- 
Riemann conditions: 



(6) 



This makes possible introducing the complex potential W = 
$ + so that 



-ic(R-R*), V = c(R + R*), R e =iR v ; 



(7) 



where R = iW* /2c, * is the complex conjugate. In terms of 
the complex function R(9, y), the dynamical and kinematic 
boundary conditions (2), (3) are as follows: 

i(R e -R* e ) + 2R e R* e +r]/c 2 = 0, y = v(0); (8) 
R(6, rfi + R*(9, r/) — 7] = B, B = const; (9) 

where B is the Bernoulli constant. It is found by averaging 
(9): B = K — 77, K = I/c is the Stokes flow, I is the wave 
impulse averaged over the period: 

I=i LJ de J y)dy = h I * {9 ' v{e)) de - 



In addition to the Laplace equation and the boundary con- 
ditions, an initial condition should be assigned. Since the 
canonical model is energy conservative, the wave total en- 
ergy can be used instead to characterise wave properties. For 
this purpose, however, the crest-to-trough height H or the 
wave steepness A = H/ A are the more convenient param- 
eters since they monotonously increase starting from linear 
waves up to the limiting configuration. Thus, using the Eqs. 
(1), (8), (9), (4) of the canonical model the following quanti- 
ties are to be found as the functions of the wave steepness A: 
the complex function R(9, y), the elevation t](9) of the free 
surface, and the wave phase speed c. 

3 The method for obtaining solutions 

3 . 1 Fractional Fourier approximations 

When working in the plane of spatial variables the solutions 
to the Laplace equation (1) in the flow domain are usually 
looked for as the following truncated Fourier series 



JV 



R(9, y) = J2^ncxp(n(y + i9)). 



(10) 



n=l 
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This approach was applied by the authors (Lukomsky et al., 
2002a,b) for calculating steep gravity waves. Fourier ex- 
pansions (10), however, become ineffective for steep waves 
with sharpening crests close to the limiting wave due to slow 
descending of the Fourier coefficients. Because of this we 
proposed (Gandzha et al., 2002; Lukomsky et al., 2002c) a 
more effective set of functions to expand the velocity poten- 
tial on the basis of the following Euler formula (see Ham- 
ming, 1962) 



n=l 



(ii) 



C™ 1 being the binomial coefficients. By choosing z(9, y) = 
exp (y — yo + iO); <r n = £„ exp (nj/o) tne following one- 
parametric expansion for the velocity potential is obtained 
after truncating the series: 



N 



R{0, y; yo) = 



Cr, 



^ (l-exp(yo-y-ifl))' 



JV 

E 



Ctr, 



i (exp(-yo) - exp (-y - i0)) n 

N 

^a„T"(0, y; y ), (12) 



T(0, y; y Q ) = (exp(-y ) - exp (-y - i6)) 



(13) 



where the normalized coefficients a n = ( n cxp(— nyo) were 
introduced to overcome infinite exponents at y — » oo. Ap- 
proximation (12) shows a formal correspondence with Pade- 
type fractional approximates. Because of this we called 
expansion (12) a "fractional Fourier expansion". It is sin- 
gular in a countable number of isolated points y = yo, 
9 = 2nk, k e Z, their location being determined by a free 
parameter y - Singular points are to be located outside the 
flow domain for calculating potential waves. At yo = oo, 
fractional Fourier expansion (12) reduces to ordinary Fourier 
expansion (10) with = (— l) n a n . Due to (11) expan- 
sions (10) and (12) are equivalent at N = oo and the con- 
vergence of (12) follows from the convergence of (10). For 
finite N and y a <~ 1, however, a fractional Fourier expansion 
converges much more rapidly than an ordinary Fourier ex- 
pansion. The reason is that a finite number of terms in (12) 
always corresponds to an infinite number of terms in (10) that 
is especially important for sharp-crested waves. 

Expansions (12), (10), and, in general, any function 
R(8, y) = R(y + i6) all satisfy the Laplace equation (1) 
exactly. The latter fact was also used by Clamond (1999, 
2003) (for finite and infinite depth, respectively) to introduce 
a renormalization principle that allows reconstructing the ve- 
locity potential in the whole domain once the velocity po- 
tential at the bottom (or any other level) is known. By ap- 
plying such renormalization to the first-order periodic solu- 
tion of KdV equation Clamond (2003) obtained the velocity 



potential being exactly the same to the first term (N = 1) 
of expansion (12), which he called a renormalized cnoidal 
wave (RCW) approximation. There may be other possibili- 
ties to improve ordinary Fourier expansion (10) besides the 
proposed fractional expansion (12). However, one should ad- 
ditionally assure the convergence of series that makes con- 
structing such generalized expansions much more difficult. 
One can see from the expansion of derivatives 

N+l 

R y (0, y; yo) = -%Re = E PnT n {9, y; y ), 

n=l 

/3 n = na n - (n - l)a„_i exp {-y ), (14) 

which follows directly from (12), that the boundary condition 
at infinite depth (4) is also satisfied exactly. 

Hereafter, only the symmetric waves are considered. In 
this case, the coefficients a n and £„ are real (in general, they 
are complex numbers for nonsymmetric waves). After taking 
into account expansions (12) and (14) the boundary condi- 
tions (8), (9) at the free surface attain the following form: 

N+l 



n i—l 

iV+1 N+l 

n i—l n 2 — 1 

N 

2 ]T a ni Re(r') 

Til- 1 



«)) =v, y = v(0); (15) 



n = B, y = rj(6). 



(16) 



Note that Eqs. (15) and (16) at N —> oo are equivalent to 
boundary conditions (8) and (9), in the class of 2tt -periodic 
functions (subharmonic waves with multiple periods are not 
taken into account in expansions (12)). 

3.2 Nonlinear transformation of the horizontal scale 

To solve boundary conditions (15) and (16) one should assign 
an appropriate approximation to the unknown elevation y = 
rj{9) of the free surface. In the plain of spatial variables, the 
Fourier series 



n(0) 



oo 

E 

n— — oo 



rjn exp(zn6»), rj^ n = rj r , 



(17) 



are often used. Taking infinite bounds is not necessary for 
profiles with rounded enough crests. However, adequate de- 
scription of sharp-crested profiles close to the limiting one re- 
quires taking into account excessively large number of modes 
due to extremely slow descending of the Fourier coefficients. 
This highly restricts practical application of (17). The fol- 
lowing non-linear transformation of the horizontal scale orig- 
inally suggested by Chen and Saffman (1980) 

#(X; 7)=X-7sinx, < 7 < 1, (18) 

allows overcoming this difficulty by stretching sharp crests to 
a more rounded configuration. As a result, the Fourier series 



M 

»7(x; 7) = E ^ ex P( in x)> v-l 



„(7) 
'In ' 



(19) 



n = 



-M 
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in the x-space with stretched crests are much more efficient 
(Lukomsky et al., 2002c). Due to nonlinear transformation 
(18) any finite number M of the coefficients ^i 7 - 1 at 7 7^ 
corresponds to infinite number of the coefficients r/ n = T]^ 
in ordinary Fourier series (7 = 0), the associated relations 
being presented in Appendix A. Thus, the role of the pa- 
rameter 7 for the series (17) in horizontal coordinate 9 is the 
same to the role of the parameter y for the series (10) in 
vertical coordinate y. 

3.3 Numerical procedure 

By means of (19) the boundary conditions (15), (16) at the 
free surface are reduced to the following system of non-linear 
algebraic equations 



V n = c 2 d n - =0, n = 0, M; (20) 

JV 

/C„ = 2^a„ 1 4" 1 )-r ? ^)=0, n = T7A; (21) 



where 



ni — l 



N+l 



N+l 



ni — 1 ri2— ni 

$m, n 2 is the Kronecker delta. The coefficients and 
ti" 1 ' " 2 ' are the Fourier harmonics of the functions Rc(T ni ) 
and Re(T ni T* ™ 2 ), respectively: 



t (n u n ^ = ±-J Rc(T ni T* " 2 ) exp(-m X )d X . (22) 


They were calculated using the fast Fourier transform (FFT). 
The Bernoulli constant is found from the kinematic equations 
(21) atn = 0: 



B=2 5> ni *<r>-„™ 



(23) 



ni — l 



The truncation of Eqs. (20), (21) was chosen for the fol- 
lowing reasons. Since the set of kinematic equations (21) 
is linear over the coefficients a n (n — 1, N), they can be 
found in terms of the harmonics rjn without using dynami- 
cal equations (20). To proceed in such a way, it is sufficient 
to take into account only the first A kinematic equations. 
Then the rest M + 1 variables c, 7t1 7) (n = 2, M) 

are found from dynamical equations (20). The last un- 
known parameter r)[ J ^ is determined by the wave steepness 
A = (??(0) — r](n)) /2ir as follows 

[(M-l)/2] 
n=l 



the square brackets designating the integer part. Since the 
wave steepness A is an integral characteristic, some wave 
properties may be missed when using it as a governing pa- 
rameter. Thus, we additionally use the first harmonic r/i of 
the elevation in the 0-space (a spectral characteristic) as an 
independent variable instead of the wave steepness. In this 
case, the first harmonic 77} in the %-space is expressed in 
terms of r)\ and the rest of the harmonics r)^ (n — 2, M) 
by means of the expression (Al) at n = 1 instead of (24). 

The set of Eqs. (20), (21) was solved by Newton's method, 
the Jacoby matrix being given in Appendix B. Starting values 
for new calculations were taken from previous runs. For large 
enough and M, the Jacoby matrix was found to become 
badly conditioned. Because of this the program realization 
was implemented in arbitrary precision computer arithmetic. 
For instance, computations at N = 150, M = 2. 5 A demand 
160-digit arithmetic that is ten times more accurate than the 
machine one. Note that such a run is equivalent in computer 
time to a run with A = 250, M = 4 A using ordinary Fourier 
approximations and because of this takes approximately 4 
times lesser computer memory. 

The truncation numbers A and M are chosen for the fol- 
lowing reasons. By fixing the number A in expansion (12) 
an approximate configuration of the velocity potential is as- 
signed. To find out a proper truncation of the series (19) 
for the elevation associated with this configuration, the num- 
ber M should be increased until the revision of solutions for 
greater M becomes less than chosen accuracy. Then the pre- 
cision to which boundary conditions (9) and (8) are satisfied 
defines the absolute errors connected with a truncation of the 
potential and elevation, respectively. The overall relative er- 
ror Er max of such an approximate solution is agreed to be the 
maximal absolute error in boundary conditions (8) and (9) all 
over the wave period divided by the Bernoulli constant B. To 
obtain a solution close to the exact one, one should gradually 
increase the number A, choosing every time a proper value 
of M, until overall desired precision is achieved. In a major- 
ity of calculations, it was sufficient to use the approximation 
M = 2.5A or lesser ones. 

The numerical scheme proposed operates with two param- 
eters yo and 7. Decreasing yo from yo = 00 to yo ~ 1 
accelerates the convergence of the fractional Fourier expan- 
sion (12) for the velocity potential, lesser A being necessary 
to retain the same accuracy. Increasing 7 from 7 = to 
7=1 — e, e^O accelerates the convergence of the expan- 
sion (19) for the elevation, lesser M being necessary to retain 
the same accuracy. These two processes, however, should 
be carried out simultaneously. Using the fractional Fourier 
expansion (12) without the transformation of the horizontal 
scale (18) was found to deteriorate the convergence of series 
(17) and, vice versa, using the transformation of the horizon- 
tal scale without the fractional Fourier expansion was found 
to deteriorate the convergence of expansion (10), with only 
slight overall benefit having been achieved. On the contrary, 
using the fractional Fourier expansion in combination with 
the nonlinear transformation of the horizontal scale proved 
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Fig. 2. The dependence of the phase speed c of steep surface waves on their steepness A. 




Fig. 3. The dependence of the phase speed c of steep surface waves on the first harmonic rji of their profile. 



Table 1. The points of extremums in phase 
of the profile for Stokes waves (N = 120, 


speed c, steepness A(e = 
M = 2.5N, y = 0.9, 7 


it A), and the first harmonic r)i 
= 0.92) 


point 


extremum 


A 


e 


c 


m 




the first max of rji 


0.1351 


0.424429 


1.0909437483 


0.1799822 


1 


the first max of c 


0.13875 


0.435896 


1.0929513818 


0.1789318 


2 


the first min of r)i 


0.14072 


0.442085 


1.0923021558 


0.1779969 


3 


the first min of c 


0.14092 


0.442713 


1.0922768392 


0.1780099 


4 


the second max of 771 


ps 0.141055 


0.443137 


as 1.092288 


0.1780222 


5* 


the second max of c 


0.141056 


0.443141 


1.0922851495 




6 


max of A (the limiting value) 


ps 0.141064 


0.443166 


ps 1.09229 


0.1780216 



'Results from Tanaka's program 
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Table 2. The wave speed c and the mean water level tjo for steep Stokes waves depending on their steepness A (e = irA)K The maximal relative errors 
Er max of corresponding approximate solutions demonstrate great advantage of fractional approximations over ordinary Fourier approximations 



Ordinary Fourier approximations 1 Fractional Fourier approximations 2 Tanaka's program 

As c t?o -Er max , % c ??o £> max , % c 



0.14 


S3 0.439823 


1.0926149034 


-2.8 


10 -20 


3.8 


10 -io 


1.0926149034° 


-2.4 


10 


-40 


2.9 ■ 


10 -20 


1.0926149034 


0.1406 


» 0.441708 


1.0923377398 


-5.1 


io- 13 


1.6 


10~ 5 


1.0923377499" 


-1.1 


10 


-22 


3.2 ■ 


io- 11 


1.0923377499 


0.14092 


« 0.442713 
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to be highly efficient (see Section 4). 
3.4 Physical quantities 

Once the coefficients a n , ry„ and the wave phase speed c are 
found, a variety of wave characteristics can be calculated. 
The velocity potential, stream function, and the horizontal 
and vertical velocities of fluid particles are as follows: 

JV 

$(6,y) = 2c5>„Im(T"(c?, y)); (25) 

71=1 



y) = 2cJ2^nRc(T n (9, y)); (26) 

71=1 

JV+1 

*e(0,v) - 2cJ2PnRc(T n (9, y)); (27) 



71=1 

JV+1 



%{6, y) = 2cJ2PnMT n (6,y)); 



(28) 



n=l 



P n = na n - (n- l)a„_i exp (-y ). 

The horizontal and vertical accelerations of fluid particles 
are as follows: 

d 2 9 

-^2 = ®08 ($0 - c) + bye^y, 

= %e ($0 - c) - $ ee $ y ; (29) 

where 

JV+2 

* oe (9,y) = -2cJ2^nlm(T n (9 7 y)); (30) 

71=1 

JV+2 

%e(9, y) = 2cJ2^nMT n (9, y)); (31) 

71=1 

Mn = n(3 n - (n- l)/3„_i exp (-y ). 



The Stokes flow i\T, the wave impulse I, and the wave ki- 
netic energy EKin are as follows (see Cokelet (1977)): 



if = 73 + 770, I = Kc, E Kin = cI/2; 



(32) 



where 770 = T]{6) is the wave mean level and is determined 
from relation (Al). 

The wave potential energy U is calculated as follows: 



u = hl\ ^ w de = l^) 2 + E W?) 2 - (33) 

" ™ 1=1 

M 

1 (<Ct) Jn) + V r>W +r7 (7) ^ 

2^0 Vi + 2-^i yVm-i+Vm+i) J- 

n-i — l 

4 Regular and irregular flows 

4. 1 Stokes waves 

The dependence c(A) of the phase speed of almost highest 
Stokes waves on their steepness calculated using fractional 
Fourier approximations is shown in Fig. 2 by the branch 1-2- 
3-6. And the corresponding dependence 0(771) of the phase 
speed on the first harmonic of the elevation is presented in 
Fig. 3 by the branch 1-2-3-4-6. The corresponding points 
of extremums in phase speed c, the first harmonic 771, and 
steepness A are presented in Table 1 . 

The advantage of fractional Fourier approximations over 
ordinary Fourier approximations is well seen from Table 2. 
There, the values of the wave phase speed c and the mean 
water level 770 calculated using these two approaches are pre- 
sented at different values of the wave steepness A up to the 
limiting one. The precise value of the mean level 770 should 
be zero for infinite depth independently of steepness. There- 
fore, the approximate values of 770 can be used to estimate the 



8 



Lukomsky and Gandzha: Fractional Fourier approximations 



=0.14092 7V=120M=300 y = °-9 7 = 0-92 

y 

0.602 
0.601 
0.6 
0.599 
0.598 
0.597 
0.596 
0.595 
0.594 

-0.006 -0.004 -0.002 0.002 0.004 0.006 

Fig. 4. The streamlines of the almost highest regular flow near the wave crest, mapped outside the domain filled by fluid, in the wave related frame of reference. 
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Fig. 5. The streamlines of the almost highest irregular flow near the wave crest, in the wave related frame of reference. The profile of the Stokes wave of the 
same steepness is also included for comparison. 
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precision of approximate results. The maximal relative er- 
rors Er max of corresponding approximate solutions are also 
presented for analysis. One can see from the relative errors 
that the benefit from using fractional Fourier approximations 
with parameters chosen varies from ten decimal orders for 
A = 0.14 (« 99.25% of the limiting steepness) to one deci- 
mal order for A — 0.141064 (almost the limiting steepness). 
After taking into account that a run using fractional Fourier 
approximations with N = 120, M = 2. 5 A needs approx- 
imately 2.5 times lesser computer time and approximately 
10 times lesser computer memory than a run using ordinary 
Fourier approximations with N = 250, M = 4A, the ad- 
vantage of fractional approximations becomes doubtless. 

Nevertheless, the results presented for the values of steep- 
ness beyond the first minimum of c (A = 0.14092) are still 
less accurate than the ones obtained from Tanaka's program, 
which are also included into Table 2 for comparison. This 
is also well seen from Fig. 2, where the second maximum 
of c (A — 0.141056, the point 5) was obtained only using 
Tanaka's program. The fractional Fourier approximation at 
yo = 0.9 is sufficient to obtain the second maximum of r\\ 
(the point 4 in Fig. 3), but not sufficient to trace the second 
maximum of c. The reason is that the value y = 0.9 used is 
optimal for the steepness corresponding to the first minimum 
of c, yet lesser y n being necessary for greater A to improve 
the precision of fractional Fourier approximations. However, 
using the present computer realization of the method the au- 
thors failed to accomplish this task due to unsatisfactory con- 
vergence of their numerical algorithm for y < 0.9. If this 
problem could be resolved fractional Fourier approximations 
in the physical plane would have a potential to become al- 
most as effective as Tanaka's method in the inverse plane. 

4.2 Irregular flows 

The flow field in Stokes waves is regular, that is, fluid par- 
ticles move slower than the wave itself all over the flow do- 
main. In the Stokes corner flow only, the fluid particle at 
the wave crest moves with velocity equal to the wave phase 
speed and, therefore, is motionless with respect to the wave. 
Because of this such points in the flow field are called the 
stagnation points. For all the Stokes waves other than the 
limiting one, a stagnation point is located outside the flow 
domain, as was at first shown by Grant (1973). The example 
of such a regular flow mapped outside the domain filled by 
fluid is presented in Fig. 4. The streamlines at the stagnation 
point O meet at right angles in accordance with the results of 
Grant (1973) and Longuet-Higgins and Fox (1978). 

Lukomsky et al. (2002b) have numerically revealed an- 
other type of waves with stagnation point located inside the 
flow domain. This stagnation point makes streamlines be dis- 
continuous in the vicinity of the wave crest and form two 
symmetric upward and downward jets. Because of this such 
waves and flows were called irregular. The example of such 
an irregular flow calculated using fractional Fourier approx- 
imations is shown in Fig. 5. Note that the streamlines at the 
stagnation point 0\ also meet at right angles. This is the gen- 



eral rule valid for all the irregular and regular flows provided 
that the stagnation point and the wave crest do not merge (see 
Appendix C for details). 

One can see from Fig. 5 that the particles from the near- 
surface layer of the irregular flow are accelerated to veloci- 
ties greater than the wave phase speed when approaching the 
crest. As a result, they form the upward jet emanating from 
the front face of the wave. The acceleration of particles at 
the base of the jet ranges from 2.5g at 9 — 0.001 to 651 at 
the wave crest. Such large accelerations of the water rising 
up the front of the wave into the jet are known to occur in 
breaking waves (see Banner and Peregrine, 1993), where the 
typical maxima in accelerations obtained from detailed un- 
steady numerical computations were reported to be around 
5g. In view of these points, irregular flows produce a sim- 
plified model for describing the initial stage of wave break- 
ing or, to be more precise, the formation of spilling breakers 
when a localized jet is formed near the crest following by 
generating whitecaps. 

Irregular wave is an approximate stationary and obviously 
unstable solution to the equations of motion. Although we 
did not investigate its unsteady evolution, the following re- 
sults are expected. The particles with velocities greater than 
the wave speed will all be jetted out away from the fluid 
and the crest will break if their is no external influx into 
the flow domain from the left downward jet. After finish- 
ing this non-stationary process the flow will become regular 
and of lesser steepness. This resembles the recurrence phe- 
nomenon observed by Longuet-Higgins and Dommermuth 
(1997) when computing unsteady non-linear development of 
the crest instabilities of almost highest Stokes waves result- 
ing in a smooth transition to a periodic wave of lower ampli- 
tude. The appearance of irregular flows may give an expla- 
nation to this phenomenon. 

Irregular flows at fixed truncation numbers N and M de- 
pend on the wave steepness A in the following way. As 
wave steepness decreases from the limiting value, the stagna- 
tion point settles down deeper into the flow domain and the 
distance between the upward and downward jets increases. 
The corresponding dependence of the phase speed of irregu- 
lar waves on their steepness c(A) is shown in Fig. 2 by the 
branch 6-8. And the dependence 0(771) of the phase speed on 
the first harmonic of the elevation is presented in Fig. 3 by the 
branch 6-7-8. While moving along the irregular branch away 
from the limiting point 6 the accuracy of approximate solu- 
tions at fixed N and M drops. Because of this the branches 
corresponding to irregular flows in Figs. 2, 3 have not yet 
stabilized with respect to increasing the truncation numbers 
N and M, although fractional approximations (12) in com- 
bination with non-linear transformation (18) are up to one 
decimal order more accurate than ordinary Fourier approx- 
imations (10) when calculating irregular waves. The loop 
in Fig. 2 still enlarges with increasing N, the cross-section 
point with the Stokes branch moving to the left. On the 
contrary, the irregular branch 6-7-8 in the dependence 0(771 ) 
(see Fig. 3) approaches to the regular branch 1-2-3-4-6 as 
accuracy is increased. They seem to coincide completely at 
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Fig. 6. The Gibbs phenomenon in the approximations to the regular Stokes wave very close to the Stokes corner flow. 
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Fig. 7. a) The behavior of the Gibbs phenomenon in the profile of the irregular wave when proceeding from the ordinary 
to fractional Fourier approximations, b) The profile of the irregular wave obtained using fractional Fourier expansions 
at different truncation numbers. The profile of the Stokes wave of the same steepness is also included for comparison. 
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Table 3. The parameters of the irregular flow at A = 0. 14092 for different approximations: c is the wave phase speed; r)i is the first harmonic of the elevation; 
rjo is the mean water level; £Y max is the maximal relative error of an approximate solution; q(0) — c is the velocity at the crest in the wave related frame of 
reference; r](0) is the height of the crest above the still water level; y s is the vertical position of the stagnation point; rj(0) — y s is the distance of the stagnation 
point from the wave crest; a m j n is the minimal slope in the Gibbs oscillations of the wave profile. The same parameters of the regular flow at A = 0.14092 
are also included for comparison* . 
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N — > 00. This may explain why irregular waves were not 
traced from stability analysis of Stokes waves. Similar to 
almost limiting Stokes waves, lesser y are necessary to im- 
prove the precision of fractional Fourier approximations and 
to stabilize the position of the irregular branch. 

Thus, our main concern is that irregular flow is an ap- 
proximate solution, whose accuracy is not sufficient enough 
to make definite conclusions even when using fractional ap- 
proximations. The only reasonable assumption is that irreg- 
ular waves correspond to singular solutions of the equations 
of motion. One can see from Fig. 5 that although streamlines 
of the irregular flow are discontinuous near the wave crest 
the wave profile rj{9) remains to be a continuous function 
everywhere. This is due to the Fourier series (19) for rj(6) 
being a single-valued smooth function, which represents an 
integral characteristic of the flow. On the contrary, the stream 
function ^(6, y) represents a local characteristic of the flow 
and is not obligatory a single-valued dependence. Because 
of this streamlines describe the discontinuous flow near the 
wave crest more adequately. It is seen from Fig. 5 that the 
profile of the irregular flow oscillates while approaching the 
wave crest, where a prominent peak (an overshoot) forms. 
This highly resembles the Gibbs phenomenon when (i) a dis- 
continuous function or (ii) a continuous function with dis- 
continuous derivatives are approximated by a truncated set 
of continuous functions (see, e.g., Arfken and Weber, 1995). 
In both cases, the Gibbs phenomenon can be used as an ex- 
cellent indicator of a singularity. 

The example corresponding to the case (i) is given in 
Appendix D, where truncated Fourier series of a function 
with infinite discontinuity are demonstrated to exhibit typi- 
cal Gibbs oscillations and an overshoot. For the case (ii), the 
limiting Stokes wave can be used as an example because of 
the sharp 120° corner at the crest corresponding to a discon- 
tinuous first derivative. Moreover, even regular Stokes waves 
very close to the Stokes corner flow also exhibit Gibbs os- 
cillations when being approximated numerically, as was re- 



vealed by Chandler and Graham (1993) from the analysis of 
Nekrasov's integral equation. The similar example obtained 
using fractional and ordinary Fourier expansions is presented 
in Fig. 6 for the Stokes wave at A = 0.14106 (w 99.9975% 
of the limiting steepness A w 0.1410635). The Gibbs phe- 
nomenon is distinctly observed even though the exact solu- 
tion is known to be not sharp-crested. The possible reason 
is that it has discontinuous derivatives higher than the first 
one being continuous. One can notice from Fig. 6 that the 
overshoot shrinks both in vertical and horizontal scales with 
increasing the accuracy of approximations in contrast to the 
example with infinite discontinuity, where the height of the 
overshoot remains almost constant (see Fig. Dl). 

Now let us analyse how an approximate irregular flow 
at fixed steepness depends on improving its numerical ac- 
curacy and make a comparison with two examples consid- 
ered above. In Fig. 7, the profiles of the irregular flow at 
A = 0.14092 are plotted for several approximations obtained 
using fractional and ordinary Fourier expansions, the corre- 
sponding flow parameters being presented in Table 3. The 
following aspects can be emphasized. 

1) The overshoot shrinks both in vertical and horizontal 
scales as the accuracy of expansions is improved, in the same 
way as for the almost limiting Stokes wave in Fig. 6. The 
horizontal and vertical dimensions of the overshoot both be- 
come approximately five times as small when using more ac- 
curate fractional approximations instead of ordinary Fourier 
approximations. And the distance r/(0) — y s between the 
wave crest and the stagnation point becomes approximately 
four times as small. This correlates with decreasing the max- 
imal relative error Er max of the solutions approximately by a 
factor of six. The conclusion can be made that the overshoot 
is likely to shrink into a single point when increasing accu- 
racy further similar to the case of the almost limiting Stokes 
wave. 

2) The overshoot has an almost fixed vertical position in 
contrast to the example with infinite discontinuity in Fig. Dl, 
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where the overshoot moves to infinity as accuracy is im- 
proved. This indicates that if the profiles of irregular flows 
are indeed discontinuous they seem to exhibit a finite discon- 
tinuity. The wave height 77(0) quickly stabilizes with increas- 
ing accuracy and seems not to tend to the height of the Stokes 
corner flow. 

3) The particle velocity at the crest q(0) — c in the wave 
related frame of reference decreases but much less rapidly 
than the relative error Er max of the corresponding solutions. 

4) The minimal slope a m ; n in the Gibbs oscillations of the 
wave profile decreases, the oscillatory tails approaching the 
horizontal line. This may be the indication that the profile 
has two symmetric finite discontinuities in the vicinity of the 
wave crest. 

These observations and the data presented do not give the 
final understanding of irregular flows at N — > 00. At present, 
we can only guess what exact solutions to the equations of 
motion irregular waves do approximate. We admit the fol- 
lowing two alternatives. 

(i) The particle velocity at the crest q(0) — c drops to zero, 
the overshoot shrinks into a single point, and the stagnation 
point merges with the wave crest. Then irregular flows ap- 
proximate a family of sharp-crested corner flows similar to 
the Stokes corner flow but of lesser steepness. 

(ii) The particle velocity at the crest q(0) — c does not tend 
to zero and the stagnation point remains in the flow domain. 
Then irregular flows are indeed discontinuous and their pro- 
files are not single-valued and have two symmetric finite dis- 
continuities in the vicinity of the wave crest. 

Further investigation is necessary to single out one of these 
possibilities or probably to come to some other conclusion. 



5 Conclusions 

Fractional Fourier approximations for the velocity potential 
in combination with non-linear transformation of the hori- 
zontal scale, which concentrates a numerical emphasis on 
the crest region, turned out to be much more efficient than 
ordinary Fourier approximations when computing both steep 
regular and irregular flows. Nevertheless, further improve- 
ment of the numerical algorithm is necessary to achieve the 
accuracy of Tanaka's method when calculating almost lim- 
iting Stokes waves and to attain the final understanding of 
irregular flows. One of the possible ways is to use the fol- 
lowing multi-term fractional expansion with several different 
parameters y^. 



K N k 



R{6, y- {?/,}) = EE 



a 



(k) 



%~o n^i ( e M-Vk) - cxp (-y - iO)) 

Although the proposed approach was formulated in the 
framework of the canonical model for infinite depth, its prac- 
tical application is much broader. Gandzha et al. (2003) suc- 
cessfully employed fractional approximations for computing 
gravity-capillary waves. When y is located inside the flow 



domain fractional approximations may be applied for calcu- 
lating vortex structures and solitary waves. The latter possi- 
bility was realized by Clamond (2003) using his renormal- 
ized cnoidal wave approximation (the first term of the frac- 
tional Fourier approximation). He computed an algebraic 
solitary wave on deep water and traced how it changes after 
taking into account surface tension (it is not known, however, 
how this algebraic solution depends on picking up higher 
terms of the fractional expansion). In the case of finite depth 
h, the fractional Fourier expansion will be as follows: 



N 

R(9, y; y Q ) = E( 



a 



(+) 



^ V (exp(-yo) - cxp (-y - h - i8)) n 

(exp(-yo) - exp (y + h + i6)) n ' 

Finally, fractional expansions may also be generalized to the 
case of 3D waves and non-ideal fluid. 

Fractional approximations allowed us to gain more de- 
tailed knowledge about the properties of irregular flows. Ir- 
regular waves were proved to correspond to singular solu- 
tions of the equations of motion. Because of this their ex- 
istence does not contradict to the uniqueness theorem of 
Garabedian (1965) since it deals with regular continuous so- 
lutions only. The following two alternatives for the exact 
solutions associated with approximate irregular flows were 
singled out: (i) the profile of an exact solution is continuous 
with sharp corner at the crest (a discontinuous first derivative) 
similar to the limiting Stokes wave but is of lesser steepness; 
(ii) the profile of an exact solution has two symmetric finite 
discontinuities in the vicinity of the wave crest and is not 
single-valued. Further analysis, however, should be carried 
out to make a final conclusion. One of the possible ways is 
to investigate how an approximate irregular flow depends on 
taking into account surface tension and to make a comparison 
with new limiting forms for gravity-capillary waves recently 
obtained by Debiane and Kharif (1996). 

To conclude note that the formation of jets from irregular 
progressive waves resembles the occurence of vertical jets 
with sharp-pointed tips from standing gravity waves forced 
beyond the maximum height, as has recently been reported 
by Longuet-Higgins (2001). 



Appendix A The relations between the Fourier coeffi- 
cients in the 9- and x-spaces 

Taking into account nonlinear transformation (18) the coeffi- 
cients in Fourier series (17) and (19) are connected as follows 

(7) (7) 
Vo = Vo -IVi , 

1 M 

Vn = -^miif^-mW-^W), (Al) 

n i — l 



n = 1, oo; J n (z) being the Bessel function of the first kind. 
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Appendix B The Jacoby matrix 

The Jacoby matrix is composed of the coefficients at the in- 
finitesimal variations 8c, Sr/^, 5r)$ {n\ = 2, M), 5a ni 
(ni = 1, AT) of the unknown variables in the following vari- 
ations of equations (20) and (21): 

N M 

SV n = a n 2 ni Sa m + E a n] ni 5l lni + 2cd n Sc; 



ni—1 
N 



ni=0 
M 



5K n = E a ili fo »i + E ^.ni^m' 



ni- 1 ni— 



where 



a 



= oi 11 +« n a 11 ^n 11 - 

n, n-i ^n—n-i 1 ^n+ni ' ^n, n ) 



JV+1 



a 



12 



ni — 1 



21 21 21 21 21 

7V+1 

= 2c 2 ^/3„ 1 (n 1 (^)-4" 1+1) e-^)- 

n i—l 

7V+1 

E (2-^ 1 ,n 2 )/3n 2 (K+^)tl"-" 2) - 
n 2 =ni 

nii (m+i lBl ) e -» -n 2 4 ni,n2+1) e- B0 )) -*n,o; 
< 2 ni = 2n 1 c 2 (4" l) -e +1) e- ! ' - 

JV+l 

2 E ^ 2 (4" 1,n2) -4" 1+1,n2) e- 



712 = 1 



Note that the variation c>?7^ should be expressed in terms of 

the rest variations 5i}n\ n = 2, M using expression (24) 
if the governing parameter is the steepness A or the relation 
(Al) at n = 1 if the governing parameter is the first harmonic 
T}\ of the elevation in the 0-space. 

Appendix C The stagnation point 

The point in the flow field, where fluid particles are mo- 
tionless in the wave related frame of reference, is called the 
stagnation point. For symmetric regular/irregular flows, the 
stagnation point is located above/below the wave crest out- 
side/inside the flow domain on the axis 9 = 0. Then its ver- 
tical position y s is determined as follows 

*o(0, y,) = c (CI) 
To find the velocity field $o(6, y), $ y (6, y) in the in- 
finitesimal vicinity 9 = 9 (9 — > 0), y = y s +y (y — > 0) of the 
stagnation point it is sufficient to linearize there the expan- 
sions (27), (28) that represent exact velocity field at N — > oo. 
For this, one should linearize the functions T n (9, y) around 
the stagnation point as follows 

T(6, y) = T(y s + y + i9) = T(y s ) + T'{y s ){y + iff); 
T n (9, y) = T n (y s ) + nT' {y s )T n -\y s ){y + iff). (C2) 



Then, after taking into account condition (CI), the Eqs. (5) 
for particle trajectories attain the following form in the vicin- 
ity of the stagnation point: 

riff rl~ °° 

^ = ay, f t =a0; a = 2c ]T n(3 n T' {y s )T n -\y s ). 

n=l 

Therefore, the equations for the streamlines are 9 = ±y. 
Actually, this is a direct consequence of the fact that any so- 
lution to the Laplace equation (1) should depend not on the 
variables (9, y) separately but on their combination y + i9. 

Thus, the streamlines meet at right angles (90°) at the stag- 
nation point. This fact is valid for any flow provided that the 
stagnation point and the wave crest do not merge. Otherwise, 
y s = 77(0); 9 and x are not independent variables, and lin- 
earization (C2) is not valid. The known example is the lim- 
iting Stokes corner, where the streamlines turn out to meet 
at 120° angle, as was at first shown by Stokes (1880). The 
question appears then what is a physical background for such 
a sharp transition from 90° to 120° angle when proceeding 
from regular Stokes flows to the Stokes corner flow? Grant 
(1973) conjectured that a continuous approach to the limit- 
ing Stokes wave is possible only if the Stokes corner flow 
has several coalescing singularities. The existence of irregu- 
lar flows may give answer where these multiple singularities 
arise from. The Stokes corner flow seems to be formed due to 
merging the stagnation points of regular and irregular flows. 



Appendix D The Gibbs phenomenon 

Consider the following 27r-periodic function 



f(x) = - In 



„ . x 
2sin- 
2 



(Dl) 



with infinite discontinuity at x = 2nk, fceZ, This func- 
tion constitutes a part of the kernel of Nekrasov's integral 
equation (see Chandler and Graham, 1993). The truncated 
Fourier series of the function f(x) have the following form 
(see Arfken and Weber, 1995): 

N 1 

fW{x) = V - cos(nx), f(x) = lim (D2) 
' — ' n jv-*oo 

n— 1 

In this case, the discontinuous function f(x) is approximated 
by the continuous functions f( N \x). One can see from 
Fig. Dl that instead of infinite discontinuities, the functions 
f( N \x) have rounded peaks (overshoots) with symmetric 
oscillatory tails that descend as the distance from the point 
of discontinuity increases. This is the well known Gibbs 
phenomenon, which always takes place when approximating 
discontinuous functions by the truncated Fourier series (see 
Arfken and Weber, 1995). As the number N is increased, 
the functions f( N >(x) approximate the function f(x) more 
precisely. The peak moves upwards and the oscillatory tails 
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Fig. Dl. The truncated Fourier series of the discontinuous function. The 
Gibbs phenomenon. 

move closer to the point of discontinuity, their amplitude and 
period decreasing. Nevertheless, the height of the peak (the 
vertical distance between the point x — and the point, 
where the oscillatory tail initiates) remains almost constant 
with increasing N. Because of this the truncated Fourier 
series representation remains unreliable in the vicinity of a 
discontinuity even for high enough N. 
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